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FOREWORD 


This report covers the work completed on the research project "Analysis 
and Computation of Internal Flow Field in a Scramjet Engine" for the period 
ending July 31, 1985. The work was supported by the NASA Langley Research 
Center (Computational Methods Branch of the High-Speed Aerodynamics Divi- 
sion) through research grant NAG-1-423; the grant was monitored by Dr. Ajay 
Kumar of the High-Speed Aerodynamics Branch. 


ABSTRACT 


The governing equations of two-dimensional chemical ly-react' 
are presented together with the global two-step chemistry model, 
nite-difference algorithm used is illustrated and the method of ci 
ing the stiffness is discussed. The computer program developed i 
solve two model problems of a premixed chemically-reacting flow, 
suits obtained are physically reasonable. 
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INVESTIGATION OF CHEMICALLY-REACTING 
SUPERSONIC INTERNAL FLOWS 


By 

T. Chitsomboon^ and S. N. Tiwari^ 

1. INTRODUCTION 

The NASA Langley Research Center has been involved in research to 
develop a supersonic combustion ramjet (scramject) propulsion system for a 
number of years (Refs. 1-3). Numerical calculations of the flow fields in 
various regions of scramjet engines have proved to be valuable in under- 
standing the nature of these flows. The flows of pure air through the 
scramjet engines were studied numerically by Kumar (Refs. 4-6). Drummond 
and Weidner (Ref. 7) assumed a complete-reaction model to study the effect 
of gaseous hydrogen fuel injection from the side walls and from the center 
strut. The complete-reaction model, however, could predict the extent of 
combustion by orders of magnitude higher than that of the actual combustion. 
A realistic combustion model of the hydrogen-air system would involve some 
60 reaction paths (Ref. 8) which would make numerical investigations very 
costly if not impossible. The global two-step combustion model (Ref. 9) 
offers an alternative between these two extremes. In thi^ global model, 
only two reaction paths and four species of reactant and oxidizer are in- 
volved. Use of this reaction model should make numerical study of the Hj- 
air reacting system feasible while maintaining important features of the 
real combustion phenomena such as the extent of mixing and the overall heat 
release. 

iGraduate Research Assistant, Department of Mechanical Engineering and 
Mechanics, Old Dominion University, Norfolk, Virginia 23508. 

^Eminent Professor, 0ep^rtment of f^chanical Engineering and Mechanics, Old 
Dominion Unive»*sity, Norfolk, Virginia 23503. 


The system of equations governing the Hj-air reacting flow is stiff due 
to the fast reaction associated with the global two-step chenistry model. 
Numerical integration of these equations by means of explicit time-dependent 
finite difference methods requires extremely small time steps to ensure 
stability. The stable time-step magnitudes are usually orders of magnitude 
smaller than the time steps dictated by the CFL number of the fluid-dynamics 
equations. The small time steps that have to be used impose a severe re- 
striction on computer resources. One way of getting around the stiff cri- 
terion is by integrating all the governing equations by f ul ly-impl icit 
methods. The system of algebraic equations resulting from ful ly-impl icit 
methods usually requires large computer storage in addition to complicated 
algorithm and long matrix-inversion time. With the advent of vector comput- 
ers, such as the VPS-32 at the NASA Langley Research Center, an explicit 
method might be more attractive because it is fully vector! zable. With the 
above mentioned stiffness, however, a ful ly-vectorized code could still 
result in a very long integration time before a steady state is reached. 

Another approach in integrating the stiff equations is by evaluating 
the terms that give rise to stiffhess implicitly while the other terms, not 
contributing to stiffness, are evaluated explicitly (Refs. 10 and 11). 
This technique allows integration of the fluid-dynamic equations by an 
explicit method using the fluid-dynamics time steps calculated from the CFL 
condition. The species equations, however, produce a block-diagonal system 
of algebraic equations since the source terms (the terms that give rise to 
the stiffness) are evaluated implicitly. The block-diagonal system can be 
inverted relatively economically as compared to inversion of a multi- 
diagonal system resulting from a ful ly-impl icit method. Also, a large 
portion of the code is still vectorizable which makes this approach 
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particularly efficient on a vector computer. Drummond et al. (Ref. 12) 
successfully used this technique, with the global two-step reaction model, 
in combination with a spectral method to solve a quasi-one-dimensional 
reacting internal flow of premixed Hj-air system. Bussing and Murman (Ref. 
13) independently used the same procedure to calculate inviscid premixed Hj- 
air reacting system by a finite-volume method. 

The purpose of this investigation is to develop a computer code to 
calculate a realistic chemically-reacting flow field in a scramjet engine by 
using the two-dimensional Navier-Stokes equations and species equations 
together with the global two-step reaction model of Ref. 8. The two-dimen- 
sional full Navier-Stokes solver of Kumar (Ref. 14-15) will be modified to 
incorporate the chemistry package. Turbulence is modeled according to the 
two-layer eddy viscosity model of Ref. 16. 
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2. ANALYSIS 


Th - overning equations of the Hj-air reacting systeiTi is presented in 
Section 2.1. Section 2.2 briefly discussed the chemistry model used in this 
study. The thermodynamics model is discussed in Section 2.3. 

2.1 Governing Equations 

The two-dimensional Navier-Stokes equations in body-fitted coordinates 
written in the strong conservation-law form can be expressed synbolically 
as 


where 
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The is determined by the Sutherland's formula whereas the is cal- 
culated by using the two-layer eddy viscosity model of Ref. 16. Pure air 
properties are used to calculate and y^ since the combustion products 
are dominated by the nitrogen species. Finally, to close the system of 
equations the ideal gas equation of state, p = pRT, is used. The gas 
constant, R, is calculated by averaging over the mixture. 

The global two-step reaction model used in this study contains five 
chemical species, i.e., Oj, HjO, Hj, OH and Nj; they will be referred to 
in the following discussion by the subscripts 1, 2, 3, 4 and 5, respec- 
t-^vely . 

The strong conservation-law form of the two-dimensional species conti- 
nuity equations can be written in body-fitted coordinates as 


3q. 3E. 3F 

+ + = JW ; i = 1,2, 3,4 (2.7) 

3t 35 3n 

where 

q. = Jpf. = Jp. (2.8a) 
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The quantities u and v are as defined in the previous discussion. The 
quantity 8 is defined as 

6 = (2.9) 

Le*Pr 


where 


Le = a/D 


( 2 . 10 ) 


and 


Pr = v/a 


(2.11) 


In this study the Lewis number (Le) is assumed to be equal to 1. This as- 
sumption is reasonable for the gaseous mixture of the present study. The 
assumption also provides a great simplification of the governing equations. 
In particula’^, the enthalpy transported by species diffusions can be com- 
bined with the heat diffusion terms to yield the terms q and q as 

X y 

indicated in Eq. (2.4). A detailed discussion can be found in Ref. 17. 
Moreover, there is no need to evaluate the species diffusion coefficient 

D. 

Note that only four species equations are actually solved. The 
nitrogen species is assumed to be inert. The mass fraction of nitrogen can 
be evaluated by the following relation; 


4 

fs = 1 - f. (2.12) 

i = l 


The above relation is nothing but the mass conservation law of the mixture. 

2.2 Chemistry Model 

The global two-step reaction model is given by 
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H 2 + 02 

^fl 

XT 

2 OH 

(2.13a) 

2 OH + H 2 

'^b2 

2 H 2 O 

(2.13b) 


where the k^'s are the forward reaction rate constants and the k^^'s are 
the backward reaction rate constants. The mmerical values of k^'s and 
kj^'s as well as the details of this chemistry model are given in Ref. 9. 

For a general reaction 
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the rate of production (or extinction) of species j, can be obtained from 
the law of mass action as 
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The rate of change of the molar concentration of species j is found by 
SLnming over a:l the reaction equations, i.e.. 
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l 
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(2.15) 


The rate of change of mass concentration of species j is found from 
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W. = C' M. 
J J J 


(2.16) 


By applying the procedure of Eqs. (2.14) through (2.16) to Eq. (2.13), the 
following mass production rates are obtained: 


^ _ '^fl‘^1^3 
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(2.17a) 

(2.17b) 

(2.17c) 

(2.17d) 


The W^- terms above are used in the right-hand side of Eq . 2.7. 

2.3 Thermodynamics Model 

The value of the specific heat at constant pressure, Cp, is assumed 
to be a linear function of temperature. The arbitrary constants are obtain- 
ed from curve fitting the temperature data of Ref. 18. For each species, 
the specific heat is, therefore, expressed as 


C . = a.T + b. 
pi 1 1 


(2.18) 


The total enthalpy of the mixture is given by 
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(2.19) 


5 T -1 

H = Z [f. / (C . dT + H.)] + i (u2 + v2) 

i^l 1 . pi 2 

Using Eq. (2.18) in Eq, (2.19) and integrating gives 
5 , 

H = Z [f.(0.5aj2 + b T + h; )1 + i (u2 + y2) (2.20) 

i=l 1 1 1 ^ 1 

The mixture gas constant can be obtained by summation over all species as 

5 

R = Z f.R. (2.21) 

i=l 

The numerical values of various thermodynamic constants are summarized in 
Table 1 below. 


Table 1. Numerical values of various thermodynamic constants. 


Species 

H*(Joule/kg) 

a 

b 

R. (Joule/kg. K) 

O2 

-271,267 

0.1198 

947 

254 

H 2 O 

-13,972,530 

0.4391 

1858 

457 


-4,200,188 

2.0546 

12867 

4035 

OH 

+1,772,591 

0.1656 

1673 

478 

f^2 

-309,483 

0.1035 

1048 

240 


The values in Table 1 above are obtained from Ref. 19. 
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3. FINITE DIFFERENCE ALGORITHM 

The Eqs. (2.1) and (2.7) can be differenced in time as (for simplicity, 
the subscript i has been dropped in the subsequent discussion) 

- q*^ = -At (— + — - + 0 (At)^ (3.1) 

3C 3n 

where (n+1) denotes the time level where solutions are being sought. It 
is to be understood that the source terms, W's, are equal to zero for the 
fluid-dynamics equations. Note that the source terms are evaluated at the 
implicit time level (n+1) to alleviate stiffness associated with fast 
chemical reactions as discussed earlier. Before developing a discrete form 
of the Eq. (3.1) the nonlinear implicit quantities must be linearized in 
order to render a linear discrete equation for solution. Using Newtonian 
linearization one obtains 

yr,+l ^ ^ (S/j (3.2) 

3q 

The q's in the above equation are the Jp^'s of the four species. Since 
botr W and q are vectors of four elenents, the term 3w/3q becomes a 
4x4 matrix, the Jacobian matrix. 

Upon defining 
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Aq" = 


3q 

n+1 n 

q -q 

3C 


- JW' 


and substituting Eq. (3.2) into Eq. (3.1) and rearranging one obtains 


(I - AtA)^ Aq^ = - At (R^) 


(3.3) 


where I is i:he 8x8 identity matrix. 

A discrete form of Eq. (3.3) can be advanced in time using tne explicit 
unsplit MacCormack algorithm (Ref. 20) as follows 


Predictor Step: 


I - AtA)" Aq'^ = - At rJ 

(3.4) 

„n+l _ , .^n 

q = q + Aq 

(3.5) 


Corrector Step: 


(I - AtA)"'^^ Aq^ = - AtRj^"'^^ 


(3.6) 


+ 1 (Aq^ + Aq^) (3.7j 
2 


The subscr'^pt f and b in Eqs. (3.4) and (3.6) denote forward and 
backward finite differences, respectively. In this study, only the steady- 
state solution is sought. The time step (At) used is determined from the 
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CFL condition: 


At = min ( ^ ^ ) (3.8) 

juj + a |v| + a 

If accurate temporal history is desired one has to choose At in a differ- 
ent manner (see, e.g., Ref. 12). 

The boundary conditions used in this study are the no-slip conditions, 
the adiabatic wall condition, the zero normal pressure gradient condition 
and the non-catalytic conditions. The non-catalytic conditions assume zero 
normal concentration gradients for all species. At the inflow boundary, all 
dependent variables are fixed according to the free stream conditions of 
interest. The method of extrapolation is used at the outflow boundary, this 
limits the application to supersonic outflow only. For a subsonic outflow, 
an appropriate method must be employed. 

Numerical smoothing functions are also added to the algorithm for both 
the fluid-dynamics and the species equations. The functions used for the 
fluid equation can be found in Refs. 4-6. For the species equations, the 
same form of the smoothing functions used in the fluid equations are 
employed; only the dependent variables are changed. The addition of the 
smoothing functions to the species equations was found to be necessary in 
order to suppress high frequency oscillations in species densities. 
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4. SOLUTION PROCEDURE 


The fluid-dynamics part of the flow is integrated in time from an ini- 
tial solution by the fully explicit MacCormack algorithm. The total 
enthalpy, H, in the energy equation has to be evaluated according to the 
formula in Eq. (2.20). The finite-difference representation of the species 
equations results in a block-diagonal system of algebraic equations and is 
inverted pointwise by an L-U decomposition procedure. 

Once all the conservative variables are obtained at a time step, the 
primitive variables have to be deciphered from them. Some difficulties 
arise in deciphering the temperature, T. The solution of the energy 
equation gives the quantity (pH-p) of the mixture. This study assumes 
that the thermodynamic property, Cp, of the mixture can he lagged one time 
step without causing any significant error. With the above assumption the 
term (pH-p) can be written as 

pH-p = oC T + — p (u2 + v^) - pRT (4.1) 

P 2 

where the starred quantity is to be evaluated at the old time level. The 
term Cp is defined as (see Eq . (2.20)) 

C = t [f, (0.5 a, T t b.)]. (4.2) 

By using Eq. (4.2) in Eq. (4.1) the temperature of the new time level can oe 
obtained directly without having to solve a second-order algebraic equation 
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in T. Finally, the pressure, p, is calculated from the ideal gas 
equation of state using the mixture gas constant R. 
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5. RESULTS AND DISCUSSION 

The non-reacting computer code of Refs. 14-15 has been modified to 
incorporate the chemistry package as described in the previous sections. 
The code is written specifically for the VPS-32 vector processing computer 
at the NASA Langley Research Center. The program was used to solve two 
model problems of supersonic combustion of premixed hydrogen-air systems. 

The geometry as well as the inflow conditions of the first test case 
are shown in Fig. 1. The equivalence rat’o (ij>) for this case is 0.2. The 
inflow gas is assumed to be a perfect mixture of hydrogen and air. A grid 
of 31x31 is used in this test case; it is shown in Fig. 2. The velocity 
vectors of the solution are illustrated in Fig. 3. The vectors are well 
behaved everywhere with indication of boundary layers developed close to the 
solid walls. Note that, in addition to the assumption of laminar flow, the 
grid use for this analysis is too coarse to resolve the detail of tne bound- 
ary layers. The main purpose of this test case was, however, to check for 
coding errors and not to capture all the details of the flow. The pressure 
contours are indicated in Fig. 4. The big bubole at the inflow boundary is 
a consequence of initial formation of the OH species which is endothermic. 
The shock wave emanating from the compression corner is clearly seen in the 
figure. Figs. 5, 6 and 7 show the concentration contours of OH, H 2 and O 2 , 
respectively. These figures exhibit one common feature in that they all 
show a concentration shock at the inflow boundary. The formation of these 
shocks are physically reasonable because the OH species was produced at a 
very fast rate. In fact, it is the fast production of the OH species that 
gives rise to the stiffness of this problem in the first place. The H 2 and 
O 2 shocks were formed simply because OH was produced at the expense of 
H 2 and O 2 . The concentration contours of H 2 O in Fig. 8, however, does 
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not indicate any concentration shock at the inflow boundary; instead HjO 
is produced in a continuous fashion. This is compatible with the physics of 
the current chemistry model in that HjO production is much slower than 
that of OH. Figs. 9 and 10 show the mass-fraction distributions of various 
species ?■' y station 15 and 5, respectively. The step changes of the 
concentration of OH, and Oj are clearly seen at the inflow boundary 
whereas continuous increase in the level of H 2 O is observed. The extent 
of reaction in Fig. 10 is higher than that of Fig. 9 due to tne fact that 
the temperature level of Fig. 10 is higher than that of Fig. 9 because tne 
flow passed through the shock wave. 

For the second test case, the geometery and inflow condition', (Fig. 11) 
are mostly the same as the first test case except is changed to 4 and 
changed to 1.0. The increase in the value of * should make this test 
case a lot more severe than the previous one. In addition, the flow is 
assumed turbulent and the grid is a 31x51 system (Fig. 12) with proper 
stretching near solid walls in order to resolve the boundary- 1 ayer's de- 

tails. '>reover, a reaction switch is incorporated into the code. Tne 
switch turns tne reaction on only when the temperature is greater than 
lOOOK. The threshold temperature of lOOOK is a realist!': one for the pres- 
ent Hj-air system. Since the inflow temperature is 900 K, there should oe 
no reaction until a shock wave raises the temperature of the mixture to 
above lOOOK. Figs. 13 and 14 show the concentration contours of OH and 
HjO, respectively. These two figures indicate, as expected, no reaction in 
the free stream; the reaction occurs only after the shock or inside the 

boundary layers where the temperature is above lOOOK. Fig. 13 indicates 

OH-concentration shock coinciding with the pressure shock whereas a more-or- 

less continuous production of HjO is observed in Fig. 14. Fig. 15 shows 
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the concentration profiles of various species along the 15th y station. 
All distributions in the figure are as expected. The OH concentration 
rises much faster than that of H 2 O whereas O 2 and H 2 concentrations 
reduce continuously. About 80 percent of is consuned within the 
resident tme of about 10x10"® second. The plot also indicates that this 
short resident time is not enough for OH to recombine with H 2 to form an 
appreciable amount of the final product, HjO. Finally, the concentration 
distributions at the second y station (inside the hot boundary layer) is 
illustrated in Fig. 16. The reaction starts right at the leading edge and 
shows a small jump near the shock wave. 



6. CONCLUSIONS 


The chemistry package added to the non-reacting Navier-Stokes code 
seems to predict the correct physics of this particular Hj-air combustion 
system. The present code assumes a premixed mixture of and air at the 
inflow boundary. For a stratified Hj^air system, a minor modification to 
the code must be made. It is observed that most of the computer time is 
used to invert the block-diagonal matrix of the species equations, A more 
efficient algorithm to invert these blocks eithe*' directly or indirectly is 
highly desirable, especially the one that is suitable for vector processing. 

Current effort is being directed toward solving a realistic supersonic 
Hj-air combustion system. Instead of a prenixed assumption, the fuel (H 2 ) 
is injected as a stream into the mainflow. Results of the current investi- 
gation will be reported at a later date. 
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Figure 1. 



Geometry and inflow conditions of test case no. 1. 
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Figure 2. Grid system of test case no. 1. 
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Figure 4 



Pressure contours of case 1. 
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Figure 5. OH-concentration contours of case 1. 


Figure 6. 



H 2 -concentrati on contours of case 1, 
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Figure /. 02-concentrat.ion contours of case 1 . 
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Figure 8. 
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HoO-concentration contours of case 1. 
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Figure 9. Mass-fraction distributions at y-station no. 15. 
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Figure 12. 31x51 grid system of case 2. 
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Figure 13. OH-concentration contours of case 2. 
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Figure 15. flass-fraction distributions 


at y-station no. 15 
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Figure 16. Mass- fraction distributions at y-station no. 15. 



